/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2009 Mathias J. Krause
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
 */

/** \file
 * Reader for Nastran 3D voxel mesh files generated by HyperMesh --
 * implementation.
 */

#include "nastranReader3D.h"

#include <iostream>
#include <math.h>
#include <fstream>
#include <sstream>
#include <cmath>

namespace olb {

/////////////// Class NastranVoxelMeshReader3D //////////////////

NastranVoxelMeshReader3D::NastranVoxelMeshReader3D(std::string& fileName,
    int offset)
  : clout(std::cout,"NastranVoxelMeshReader3D") {
  _fileName = fileName;
  _offset = offset;

  // Vector containing the lines of file _fileName
  std::vector<std::string> line_v = linesOfFile2Vector(_fileName);
  std::vector<std::string>::iterator line_iter = line_v.begin();

  // A pointID is mapped to its coordnates
  std::map<int, std::vector<double> > pointID2data;
  // A quadID is mapped to its pointIDs
  std::map<int, std::vector<int> > quadID2data;
  // A hexaID is mapped to its pointIDs and colour
  std::map<int, std::vector<int> > hexaID2data;

  while (line_iter != line_v.end()) {

    //Map the grid points
    writeGridPoint(line_iter, pointID2data);

    //Map the quads
    //writeQuad(line_iter, quadID2data );

    //Map the hexas
    writeHexa(line_iter, hexaID2data);
    //Get color
    writeHexaColor(line_iter, hexaID2data);

    line_iter++;
  }

  computeMotherCuboid(pointID2data, hexaID2data);

  computeVoxelColour(pointID2data, hexaID2data);
}

NastranVoxelMeshReader3D::~NastranVoxelMeshReader3D() {

}

void NastranVoxelMeshReader3D::writeVtiFile() {

  /// Matrix of voxels with material number / boundary number
  olb::ScalarField3D<unsigned short>* geometryData;

  geometryData = new olb::ScalarField3D<unsigned short>(_nx, _ny, _nz);
  geometryData->construct();

  for (int i = 0; i < _nx; i++) {
    for (int j = 0; j < _ny; j++) {
      for (int k = 0; k < _nz; k++) {
        geometryData->get(i, j, k) = 0;
      }
    }
  }

  std::map<std::vector<unsigned short>, unsigned short>::iterator iter2;
  for (iter2 = _voxel2colour.begin(); iter2 != _voxel2colour.end(); iter2++) {
    geometryData->get((iter2->first)[0], (iter2->first)[1],
                      (iter2->first)[2]) = iter2->second;
  }

  std::string fileNameVti = _fileName;
  olb::VtkImageOutput3D<unsigned short> vtkOut(fileNameVti, 1);
  vtkOut.writeData<unsigned short,unsigned short> (*geometryData, "geometry");
  delete geometryData;
}

std::vector<std::string> NastranVoxelMeshReader3D::linesOfFile2Vector(
  std::string fileName) {

  std::ifstream infile;
  std::string line;
  std::vector<std::string> line_v;

  infile.open(fileName.c_str());
  if (!infile)
    clout << "Error opening input file" << std::endl;

  while (getline(infile, line)) {
    line_v.push_back(line);
  }
  return line_v;
}

template<>
double NastranVoxelMeshReader3D::extractValue<double>(std::string str,
    int begin, int length) {

  std::string word;
  word = str.substr(begin, length);
  double value;
  std::stringstream sstr;
  sstr << word;
  sstr >> value;
  return value;
}

template<>
int NastranVoxelMeshReader3D::extractValue<int>(std::string str, int begin,
    int length) {

  std::string word;
  word = str.substr(begin, length);
  int value;
  std::stringstream sstr;
  sstr << word;
  sstr >> value;
  return value;
}

void NastranVoxelMeshReader3D::writeGridPoint(
  std::vector<std::string>::iterator &lineIter, std::map<int,
  std::vector<double> > &pointID2data) {

  std::string new_string(*lineIter);
  if (new_string.find("GRID*", 0) != std::string::npos) {

    int gridID;
    std::vector<double> gridCo;

    // Get grid ID
    gridID = extractValue<int> (new_string, 8, 16);

    // Get grid x coordinate
    gridCo.push_back(extractValue<double> (new_string, 40, 16));

    // Get grid y coordinate
    gridCo.push_back(extractValue<double> (new_string, 56, 16));

    lineIter++;
    new_string = *lineIter;

    // Get grid z coordinate
    gridCo.push_back(extractValue<double> (new_string, 8, 16));

    pointID2data[gridID] = gridCo;
  }
}

void NastranVoxelMeshReader3D::writeQuad(
  std::vector<std::string>::iterator &lineIter, std::map<int,
  std::vector<int> > &quadID2data) {

  std::string new_string(*lineIter);
  if (new_string.find("CQUAD4", 0) != std::string::npos && new_string.find(
        "$", 0) == std::string::npos) {

    int quadID;
    std::vector<int> quadCo;

    //Get quad ID
    quadID = extractValue<int> (new_string, 8, 8);

    //Get quad points 1 - 4
    quadCo.push_back(extractValue<int> (new_string, 16, 8));
    quadCo.push_back(extractValue<int> (new_string, 24, 8));
    quadCo.push_back(extractValue<int> (new_string, 32, 8));
    quadCo.push_back(extractValue<int> (new_string, 40, 8));

    quadID2data[quadID] = quadCo;
  }
}

void NastranVoxelMeshReader3D::writeHexa(
  std::vector<std::string>::iterator &lineIter, std::map<int,
  std::vector<int> > &hexaID2data) {

  std::string new_string(*lineIter);
  if (new_string.find("CHEXA", 0) != std::string::npos && new_string.find(
        "$", 0) == std::string::npos) {

    int hexaID;
    std::vector<int> hexaCo;

    //Get hexa ID
    hexaID = extractValue<int> (new_string, 8, 8);

    //Get hexa points 1 - 6
    hexaCo.push_back(extractValue<int> (new_string, 24, 8));
    hexaCo.push_back(extractValue<int> (new_string, 32, 8));
    hexaCo.push_back(extractValue<int> (new_string, 40, 8));
    hexaCo.push_back(extractValue<int> (new_string, 48, 8));
    hexaCo.push_back(extractValue<int> (new_string, 56, 8));
    hexaCo.push_back(extractValue<int> (new_string, 64, 8));

    lineIter++;
    new_string = *lineIter;

    //Get hexa points 7 + 8
    hexaCo.push_back(extractValue<int> (new_string, 8, 8));
    hexaCo.push_back(extractValue<int> (new_string, 16, 8));

    //initialize the color (with zero)
    hexaCo.push_back(0);

    hexaID2data[hexaID] = hexaCo;
  }
}

void NastranVoxelMeshReader3D::writeHexaColor(
  std::vector<std::string>::iterator &lineIter, std::map<int,
  std::vector<int> > &hexaID2data) {

  std::string new_string(*lineIter);

  if (new_string.find("$HMMOVE", 0) != std::string::npos) {

    int voxelIDb4thru = -1;
    int color = extractValue<int> (new_string, 8, 8);
    int temp = 0;

    lineIter++;
    new_string = *lineIter;

    while (new_string.length() > 8) {
      int ll = new_string.length() - 1;
      for (int i = 8; i < ll; i += 8) {
        if (new_string.substr(i, 8) != "THRU    ") {
          temp = extractValue<int> (new_string, i, 8);
          hexaID2data[temp][8] = color;
          if (voxelIDb4thru != -1) {
            for (int j = voxelIDb4thru + 1; j < temp; j++) {
              hexaID2data[j][8] = color;
            }
            voxelIDb4thru = -1;
          }
        } else {
          voxelIDb4thru = temp;
        }
      }
      lineIter++;
      new_string = *lineIter;
    }
  }
}

void NastranVoxelMeshReader3D::computeMotherCuboid(std::map<int, std::vector<
    double> > &pointID2data, std::map<int, std::vector<int> > &hexaID2data) {

  double xM, yM, zM;
  int counter = 1;
  _h = 0;

  while (!(_h > 0)) {
    if (pointID2data[(hexaID2data.begin()->second)[0]][0]
        == pointID2data[(hexaID2data.begin()->second)[counter]][0]
        && pointID2data[(hexaID2data.begin()->second)[0]][1]
        == pointID2data[(hexaID2data.begin()->second)[counter]][1]) {
      _h = pointID2data[(hexaID2data.begin()->second)[0]][2]
           - pointID2data[(hexaID2data.begin()->second)[counter]][2];
      if (_h < 0)
        _h = -_h;
    }
    counter++;
  }

  _x0 = (pointID2data.begin()->second)[0];
  xM = (pointID2data.begin()->second)[0];
  _y0 = (pointID2data.begin()->second)[1];
  yM = (pointID2data.begin()->second)[1];
  _z0 = (pointID2data.begin()->second)[2];
  zM = (pointID2data.begin()->second)[2];

  std::map<int, std::vector<double> >::iterator iter;
  for (iter = pointID2data.begin(); iter != pointID2data.end(); iter++) {
    if ((iter->second)[0] < _x0)
      _x0 = (iter->second)[0];
    if ((iter->second)[0] > xM)
      xM = (iter->second)[0];
    if ((iter->second)[1] < _y0)
      _y0 = (iter->second)[1];
    if ((iter->second)[1] > yM)
      yM = (iter->second)[1];
    if ((iter->second)[2] < _z0)
      _z0 = (iter->second)[2];
    if ((iter->second)[2] > zM)
      zM = (iter->second)[2];
  }

  //rounding correctly
  _nx = int((xM - _x0) / _h + 2 * _offset + 0.5);
  _ny = int((yM - _y0) / _h + 2 * _offset + 0.5);
  _nz = int((zM - _z0) / _h + 2 * _offset + 0.5);

  //_x0-=_offset*_h;
  //_y0-=_offset*_h;
  //_z0-=_offset*_h;
}

void NastranVoxelMeshReader3D::computeVoxelColour(std::map<int, std::vector<
    double> > &pointID2data, std::map<int, std::vector<int> > &hexaID2data) {

  std::map<int, std::vector<int> >::iterator iter2;
  for (iter2 = hexaID2data.begin(); iter2 != hexaID2data.end(); iter2++) {
    double ix = (pointID2data[(iter2->second)[0]])[0];
    double iy = (pointID2data[(iter2->second)[0]])[1];
    double iz = (pointID2data[(iter2->second)[0]])[2];

    for (int i = 1; i < 8; i++) {
      if (ix > (pointID2data[(iter2->second)[i]])[0])
        ix = (pointID2data[(iter2->second)[i]])[0];
      if (iy > (pointID2data[(iter2->second)[i]])[1])
        iy = (pointID2data[(iter2->second)[i]])[1];
      if (iz > (pointID2data[(iter2->second)[i]])[2])
        iz = (pointID2data[(iter2->second)[i]])[2];
    }
    //rounding correctly
    std::vector<unsigned short> voxelCo;
    voxelCo.push_back(int((ix - _x0) / _h + 0.5) + _offset);
    voxelCo.push_back(int((iy - _y0) / _h + 0.5) + _offset);
    voxelCo.push_back(int((iz - _z0) / _h + 0.5) + _offset);

    _voxel2colour[voxelCo] = (iter2->second)[8];
  }
}

void NastranVoxelMeshReader3D::setColour(unsigned short iX, unsigned short iY,
    unsigned short iZ, unsigned short colour) {

  std::vector<unsigned short> co;
  co.push_back(iX);
  co.push_back(iY);
  co.push_back(iZ);

  if (colour == 0) {
    _voxel2colour.erase(_voxel2colour.find(co));
    return;
  } else {
    _voxel2colour[co] = colour;
    return;
  }
}

unsigned short NastranVoxelMeshReader3D::getColour(unsigned short iX,
    unsigned short iY, unsigned short iZ) {

  std::vector<unsigned short> co;
  co.push_back(iX);
  co.push_back(iY);
  co.push_back(iZ);

  std::map<std::vector<unsigned short>, unsigned short>::iterator foundIter;
  foundIter = _voxel2colour.find(co);

  if (foundIter == _voxel2colour.end())
    return 0;
  else
    return foundIter->second;
}

void NastranVoxelMeshReader3D::addBoundaryLayer(unsigned short from,
    unsigned short to, unsigned short closeTo) {

  for (int i = _offset; i < _nx - _offset; i++) {
    for (int j = _offset; j < _ny - _offset; j++) {
      for (int k = _offset; k < _nz - _offset; k++) {
        if (getColour(i, j, k) == closeTo) {
          if (getColour(i - 1, j - 1, k + 1) == from)
            setColour(i - 1, j - 1, k + 1, to);
          if (getColour(i - 1, j + 1, k - 1) == from)
            setColour(i - 1, j + 1, k - 1, to);
          if (getColour(i + 1, j - 1, k - 1) == from)
            setColour(i + 1, j - 1, k - 1, to);
          if (getColour(i + 1, j + 1, k - 1) == from)
            setColour(i + 1, j + 1, k - 1, to);
          if (getColour(i - 1, j + 1, k + 1) == from)
            setColour(i - 1, j + 1, k + 1, to);
          if (getColour(i + 1, j - 1, k + 1) == from)
            setColour(i + 1, j - 1, k + 1, to);
          if (getColour(i - 1, j - 1, k - 1) == from)
            setColour(i - 1, j - 1, k - 1, to);
          if (getColour(i + 1, j + 1, k + 1) == from)
            setColour(i + 1, j + 1, k + 1, to);

          if (getColour(i + 1, j - 1, k) == from)
            setColour(i + 1, j - 1, k, to);
          if (getColour(i - 1, j - 1, k) == from)
            setColour(i - 1, j - 1, k, to);
          if (getColour(i + 1, j + 1, k) == from)
            setColour(i + 1, j + 1, k, to);
          if (getColour(i - 1, j + 1, k) == from)
            setColour(i - 1, j + 1, k, to);
          if (getColour(i - 1, j, k - 1) == from)
            setColour(i - 1, j, k - 1, to);
          if (getColour(i + 1, j, k - 1) == from)
            setColour(i + 1, j, k - 1, to);
          if (getColour(i - 1, j, k + 1) == from)
            setColour(i - 1, j, k + 1, to);
          if (getColour(i + 1, j, k + 1) == from)
            setColour(i + 1, j, k + 1, to);
          if (getColour(i, j - 1, k - 1) == from)
            setColour(i, j - 1, k - 1, to);
          if (getColour(i, j + 1, k - 1) == from)
            setColour(i, j + 1, k - 1, to);
          if (getColour(i, j - 1, k + 1) == from)
            setColour(i, j - 1, k + 1, to);
          if (getColour(i, j + 1, k + 1) == from)
            setColour(i, j + 1, k + 1, to);

          if (getColour(i - 1, j, k) == from)
            setColour(i - 1, j, k, to);
          if (getColour(i + 1, j, k) == from)
            setColour(i + 1, j, k, to);
          if (getColour(i, j - 1, k) == from)
            setColour(i, j - 1, k, to);
          if (getColour(i, j + 1, k) == from)
            setColour(i, j + 1, k, to);
          if (getColour(i, j, k - 1) == from)
            setColour(i, j, k - 1, to);
          if (getColour(i, j, k + 1) == from)
            setColour(i, j, k + 1, to);
        }
      }
    }
  }
}

//template <>
//void NastranVoxelMeshReader3D::readNastranFile(/*CuboidGeometry3D<T>* _cuboidGeometry*/) {
//
//    //CuboidGeometry3D<T> cuboidGeometry(_x0, _y0, _z0, _h, _nx, _ny, _nz, 10000);
//    CuboidGeometry3D<double> cuboidGeometry(0, 0, 0, 1, _nx, _ny, _nz, 10000);
//
//    cuboidGeometry.printStatistics();
//
//    int nC = cuboidGeometry.get_nC();
//    std::vector<Cuboid3D<double> > cuboids;
//    std::vector<bool> allZero;
//
//    for (int iC = 0; iC < nC + 1; iC++) {
//        allZero.push_back(true);
//    }
//
//    std::map<std::vector<unsigned short>, unsigned short>::iterator iter2;
//    for (iter2 = _voxel2colour.begin(); iter2 != _voxel2colour.end(); iter2++) {
//        if (iter2->second != 0) {
//            allZero[cuboidGeometry.get_iC((iter2->first)[0], (iter2->first)[1],
//                    (iter2->first)[2])] = false;
//        }
//    }
//
//    if (!allZero[nC])
//        clout
//                << "ERROR: There is at least one node outside the CuboidGeometry!"
//                << std::endl;
//
//    for (int iC = 0; iC < nC; iC++) {
//        if (allZero[iC])
//            cuboidGeometry.remove(iC);
//    }
//
//    cuboidGeometry.printStatistics();
//}

void NastranVoxelMeshReader3D::readNastranFile(BlockGeometry3D* blockGeometry) {

  /// Matrix of voxels with material number / boundary number
  olb::ScalarField3D<unsigned short>* geometryData;

  geometryData = new olb::ScalarField3D<unsigned short>(_nx, _ny, _nz);
  geometryData->construct();

  for (int i = 0; i < _nx; i++) {
    for (int j = 0; j < _ny; j++) {
      for (int k = 0; k < _nz; k++) {
        geometryData->get(i, j, k) = 0;
      }
    }
  }

  std::map<std::vector<unsigned short>, unsigned short>::iterator iter2;
  for (iter2 = _voxel2colour.begin(); iter2 != _voxel2colour.end(); iter2++) {
    geometryData->get((iter2->first)[0], (iter2->first)[1],
                      (iter2->first)[2]) = iter2->second;
  }

  blockGeometry->reInit(_x0 + _offset * _h, _y0 + _offset * _h, _z0 + _offset
                        * _h, _h, _nx - 2 * _offset, _ny - 2 * _offset, _nz - 2 * _offset,
                        _offset, geometryData);
  delete geometryData;
}

} // namespace olb
